#include<iostream>
#include<limits>
#include<cmath>

class Function
{
public:
  double operator()(double _x);
  double diff(double _x);
};

class EquationSolver
{
public:
  virtual double solve() = 0;
};

class Bisection : public EquationSolver
{
private:
  double a;
  double b;
  double delta;
  int M;
  Function f;
public:
  Bisection(double _a,
	    double _b,
	    double _delta,
	    int _M,
	    Function _f);
  double solve();
};

Bisection::Bisection(double _a, double _b, double _delta, int _M, Function _f)
{
  a = _a;
  b = _b;
  delta = _delta;
  M = _M;
  f = _f;
}

double Bisection::solve()
{
  double u = f(a);
  double v = f(b);
  int k;
  double h;
  double c;
  double w;
  for (k = 1;k <= M; k++)
    {
      h = b - a;
      c = a + h/2;
       w = f(c);
      if( fabs(h) < delta || fabs(w) < std::numeric_limits<double>::epsilon() )
	{
	  break;
	}
      else if(w *u < 0)
	{
	  b = c;
	  v = w;
	}
      else
	{
	  a = c;
	  u = w;
	}
    }
  std::cout << "c="<< c << std::endl;
  std::cout << "h="<< h << std::endl;
  std::cout << "k="<< k << std::endl;
  return 0;
}

class Newton : public EquationSolver
{
private:
  double x0;
  int M;
  Function f;
public:
  Newton(
	 double _x0,
	 int _M,
	 Function _f);
  double solve();
};

Newton::Newton(double _x0, int _M, Function _f)
{
  x0 = _x0;
  M = _M;
  f = _f;
}

double Newton::solve()
{
  double x = x0;
  int k;
  for (k = 0; k <= M; k++ )
    {
      double u = f(x);
      if ( fabs(u) < std::numeric_limits<double>::epsilon())
	{
	  break;
	}
      x = x - u/f.diff(x);
    }
  std::cout << "x="<< x << std::endl;
  std::cout << "k="<< k << std::endl;
  return 0;
}


class Secant : public EquationSolver
{
private:
  double x0;
  double x1;
  double delta;
  int M;
  Function f;
public:
  Secant(
	 double _x0,
	 double _x1,
	 double _delta,
	 int _M,
	 Function _f);
  double solve();
};

Secant::Secant(double _x0,double _x1,double _delta,int _M, Function _f)
{
  x0 = _x0;
  x1 = _x1;
  delta = _delta;
  M = _M;
  f = _f;
}

double Secant::solve()
{
  double x_n = x1;
  double x_n1 = x0;
  double x_n2;
  double w;
  int k;
  double u = f(x_n);
  double v = f(x_n1);
  for (k = 2; k <= M; k++ )
    {
      if (fabs(u) > fabs(v))
	{
	  x_n2 = x_n1;
	  x_n1 = x_n;
	  x_n = x_n2;
	  w = u;
	  u = v;
	  v = w;
	}
      double s = (x_n - x_n1) / (u - v);
      x_n1 = x_n;
      v = u;
      x_n = x_n - u * s;
      u = f(x_n);
      if (fabs(x_n - x_n1) < delta || fabs(u) < std::numeric_limits<double>::epsilon())
	{
	  break;
	}
    }
  std::cout << "x_n="<< x_n << std::endl;
  std::cout << "x_n1="<< x_n1 << std::endl;
  std::cout << "k="<< k << std::endl;
  return 0;
}

